function [istp,cc,h,nnx]=phase_plot(istp)
% sumstep has been re-written as a function. 
% input = time step you want to watch
%
%
%istp=20;
xmin = -30;
xmax = 30;

zmin = -20;
zmax = 4;

rho = 3300;
alpha = 3e-5;

if exist('platform.s','file') == 0
  plat = 'native';
elseif exist('platform.s','file') == 2
  plat = 'ieee-le';
end
fid_t = fopen('time.0','r',plat);
fid1 = fopen('sxx.0','r',plat);
fid2 = fopen('syy.0','r',plat);
fid3 = fopen('srII.0','r',plat);
fid4 = fopen('temp.0','r',plat);
fid5 = fopen('phas.0','r',plat);

loadflac;
tsp = round(t(istp)/3600/24/365/1000);
nnx=nx;
for ii = 1:length(t)
    for jj = 1:nx
        %dZ = abs(diff([reshape(meshzz(1:end-1,jj,ii),ny,1);meshzz(1,jj,ii)+meshzz(end,1,1)]));
        dZ = -(diff([reshape(meshzz(1:end-1,jj,ii),ny,1);meshzz(1,jj,ii)+meshzz(end,1,1)]));
        drho = reshape(rho*(1-(alpha)*A4(:,jj,ii)),1,ny);
        tmass(jj) = drho*dZ*1e3;
    end
    zbox = meshzz(1,1,1)-meshzz(end,1,1);
    if tmass==0
        ziso(ii,:) = zeros(size(tmass));
    else
        ziso(ii,:) = zbox./(tmass./tmass(1))-zbox;
    end
end

%figure
clf
kt=istp;
a1=zeros(ny,nx); 

if istp > 3
    ma1 = (A1(:,:,kt-3) + A1(:,:,kt-2) + A1(:,:,kt-1) + A1(:,:,kt))./4;
    ma3 = (A3(:,:,kt-3) + A3(:,:,kt-2) + A3(:,:,kt-1) + A3(:,:,kt))./4;
    ma4 = (A4(:,:,kt-3) + A4(:,:,kt-2) + A4(:,:,kt-1) + A4(:,:,kt))./4;
else
    ma1=A1(:,:,kt);ma3=A3(:,:,kt);ma4=A4(:,:,kt);
end
ma5 = A5(:,:,kt);
ma5=round(ma5); % rounding the phase numbers
xcnr=zeros(ny*nx,4); zcnr=xcnr;
xcnr(:)=Xcnr(:,:,kt);  zcnr(:)=Zcnr(:,:,kt);

vfac=4;qsize=1;
jv=[1:vfac:nx]; iv=[1:vfac:ny];

Vx=zeros([length(iv) length(jv)]); Vy=Vx; Xv=Vx; Yv=Vx;
Vx(:)=vx(iv,jv,kt); Vy(:)=vz(iv,jv,kt);
Xv(:)=meshxx(iv,jv,kt); Yv(:)=meshzz(iv,jv,kt);





 [c,h]=contour(Xcnt(:,:,kt),Zcnt(:,:,kt),A4(:,:,kt),[600 600],'k','linewidth',2);
 
 gg=length(c(1,:));
 cc(1,:)=c(1,2:gg);
 cc(2,:)=c(2,2:gg);
% contour(Xcnt(:,:,kt),Zcnt(:,:,kt),A4(:,:,kt),[600 600],'k')%,'linewidth',2)

figure;
patch(xcnr',zcnr',reshape(ma5',length(xcnr),1)','edgecolor','none')
 hold on;
 plot(cc(1,:),cc(2,:),'k','linewidth',2)

%hold on; contour(Xp,Zp,a4,[600 600],'k','linewidth',2)
%plot([-100 100], [-6 -6],'color',[.8 .8 .8],'linewidth',2)
%set(gca,'layer','top','box','on')
%caxis([0 1300])
%[pos0,pos]=fixposition(1,0,1,0.05,.5,0,1,0);
  cb=colorbar;
  cbxlb = get(cb,'Xlabel');
%cbpos=get(cb,'Position');
%set(cb,'Position',[cbpos(1)+.12 cbpos(2) cbpos(3)*.5 cbpos(4)*.8])
  %set(cbxlb,'string','Phase')
xlabel('Across-Axis Distance (km)')
ylabel('Depth (km)')
axis([xmin xmax zmin zmax])
axis image

title(['Time = ',num2str(t(kt)/3600/24/365/1000),' kyrs'])
%axis([xmin xmax -25 5])
%set(gca,'XTick',-45:15:45)
%set(gca,'Position',pos)
